Loop algorithm for classical Heisenberg models with spin-ice type degeneracy 
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In many frustrated Ising models, a single-spin flip dynamics is frozen out at low temperatures com- 
pared to the dominant interaction energy scale because of the discrete "multiple valley" structure of 
degenerate ground-state manifold. This makes it difficult to study low-temperature physics of these 
frustrated systems by using Monte Carlo simulation with the standard single-spin flip algorithm. A 
typical example is the so-called spin ice model, frustrated ferromagnets on the pyrochlore lattice. 
The difficulty can be avoided by a global-flip algorithm, the loop algorithm, that enables to sample 
over the entire discrete manifold and to investigate low-temperature properties. We extend the loop 
algorithm to Heisenberg spin systems with strong easy-axis anisotropy in which the ground-state 
manifold is continuous but still retains the spin-ice type degeneracy. We examine different ways of 
loop flips and compare their efficiency. The extended loop algorithm is applied to the following two 
models, a Heisenberg antiferromagnet with easy-axis anisotropy along the z axis, and a Heisenberg 
spin ice model with the local (111) easy-axis anisotropy. For both models, we demonstrate high 
efficiency of our loop algorithm by revealing the low-temperature properties which were hard to 
access by the standard single-spin flip algorithm. For the former model, we examine the possibility 
of order-from-disorder and critically check its absence. For the latter model, we elucidate a gas- 
liquid-solid transition, namely, crossover or phase transition among paramagnet, spin-ice liquid, and 
ferromagnetically-ordered ice-rule state. 

PACS numbers: 75.10.-b, 75.10.Hk, 75.40.Mg 



I. INTRODUCTION 

Geometrically frustrated systems have attracted much 
attention because of fascinating phenomena arising from 
competing interactions^. Frustration prevents simul- 
taneous optimization of all interaction energies, which 
suppresses long-range ordering and may lead to a 
macroscopic number of energetically (nearly-)degenerate 
ground states. Such ground-state manifold plays a deci- 
sive role in low-temperature(T) physics under the influ- 
ence of quantum/thermal fluctuations and external per- 
turbations. It is highly important to clarify the structure 
of the manifold and to take the statistical average over 
the entire manifold for understanding low-T properties 
in frustrated systems. 

As a typical example, we consider antiferromagnets 
with classical spins on the pyrochlore lattice. The py- 
rochlore lattice is a three-dimensional frustrated struc- 
ture given by a corner-sharing network of tetrahedra, as 
shown in Fig. [TJ When the system has the Heisenberg 
0(3) symmetry and the exchange interaction is limited 
to nearest-neighbor sites, any long-range ordering does 
not occur and the ground-state manifold has continuous 
macroscopic degeneracy^. The manifold is identified by a 
collection of local constraints, that is, the summation of 
spin vectors on four vertices should vanish in every tetra- 
hedron. This condition is underconstraint and leaves two 
angles undetermined in each tetrahedron, resulting in the 
continuous macroscopic degeneracy^. A similar zero- 
sum local constraint can be found in an Ising ferromagnet 



on the pyrochlore lattice with local cubic (111) axes, the 
so-called spin ice model"*, which is equivalent to a py- 
rochlore Ising antiferromagnet with a global anisotropy 
axis 7 -. In this case, the local constraint enforces two spins 
pointing inward and two spins pointing outward in every 
tetrahedron, as exemplified in Fig. [I] This two- in two-out 
constraint is called the ice rule because of an analogy to 
the constraint on positions of protons in hexagonal ice&^. 
The ice rule is also underconstraint, leading to the disor- 
dered ground-state with macroscopic degeneracy. In this 
case, the manifold has a discrete nature because of the 
Ising spin degree of freedom. Intermediate type mani- 
folds, namely not discrete but not fully 0(3), also ap- 
pear in variants of pyrochlore antiferromagnets, such as 
anisotropic Heisenberg models^ and bilinear-biquadratic 
Heisenberg models^. 

In the nearest-neighbor antiferromagnetic Heisenberg 
model, all energetically-degenerate spin configurations 
in the continuous manifold are connected by continuous 
changes of spin directions without energy cost because 
of the continuous structure of the manifold. This is in- 
deed observed in classical Monte Carlo (MC) studies of 
the pyrochlore Heisenberg antiferromagnets; a standard 
single-spin flip update is efficient to sample over the en- 
tire manifold down to very low T compared to the ex- 
change energy scale J. In contrast, for discrete Ising-type 
or continuous but strongly anisotropic manifolds, degen- 
erate spin configurations are separated by large energy 
barriers of the order of J, and a single-spin flip does not 
work at low T <C J ■ In fact, in the spin ice model with 
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FIG. 1. (color online). The pyrochlore lattice composed of 
a three-dimensional network of corner-sharing tetrahedra. A 
16-site cubic unit cell is shown. A and B represent upward and 
downward tetrahedra, respectively. Ising spins are denoted 
by arrows. Each spin axis is along the local (111) axis, which 
goes from the site to the center of a neighboring tetrahedron. 
Black (filled) circles represent spins pointing inward, while 
white (open) circles represent spins pointing outward in terms 
of type-A tetrahedra. The hexagon with a bold dashed line 
denotes one of the shortest loops on which a flip of all spins 
(colors) transforms an ice-rule state to another ice-rule state. 
See the text for details. 



long-range dipole interactions, it is hard to clarify low-T 
properties by single-spin flip MC calculations, and there 
was a controversy about the possibility of long-range or- 
dering^—. 

To overcome such difficulty coming from "multiple val- 
ley" structure of the degenerate manifold, it is neces- 
sary to consider a global flip which connects different de- 
generate configurations. For the spin ice models, this 
is achieved by introducing the loop-flip algorithm, in 
which one reverses all Ising spins on a specific closed 
loop passing through tetrahedra 1 ^— ; the loop is chosen 
so that the Ising spins are inward and outward alterna- 
tively along the loop (in/outward is defined on one of 
two different types of tetrahedra, say, type-A tetrahedra 
in Fig. [IJ. This loop flip enables to transform an ice- 
rule state to another ice-rule state bypassing the energy 
barriers. Melko et al. applied the loop algorithm and suc- 
cessfully observed symmetry breaking emergent from the 
ice-rule manifold in the dipolar spin-ice model without 
severe dynamical freezin g 14 ' 17 . The loop algorithm has 
been successfully applied to study the low-T properties 
of spin-ice type Ising models 1 ^— . 

The loop algorithm is well defined for systems with dis- 
crete Ising-type spins. However, when the discreteness is 
relaxed and spins can fluctuate around the anisotropy 
axis, it becomes nontrivial how to define the loop with 
alternating in/outward spins. Moreover, it is also un- 
clear how thermal spin fluctuations affect the acceptance 



rate of the loop flip update. As long as the degenerate 
manifold retains a "multiple valley" structure with large 
energy barriers, a single-spin flip update becomes ineffi- 
cient and some global flip is indispensable for taking the 
statistical average in an ergodic way. Such problems are 
encountered in many systems, e.g., frustrated Heisenberg 
models with strong easy-axis anisotropy or with large bi- 
quadratic interactions. To elucidate low-T properties in 
this type of frustrated systems, it is desired to establish a 
global flip update applicable to systems with continuous 
but ice-rule type degenerate manifold. Note that the sit- 
uation has an aspect similar to the Wolff's extension^! of 
the Swendsen-Wang cluster algorithm^ for conventional 
unfrustrated magnets. 

In this paper, we extend the loop algorithm to Heisen- 
berg spin systems with strong easy-axis anisotropy in 
which the ground-state manifold is continuous but retains 
the ice-rule type "multiple valley" structure. We describe 
how to define the loop in such systems, and discuss two 
different types of loop flip. Because of the continuous de- 
grees of freedom of Heisenberg spins, the way of flipping 
is not unique and its efficiency depends on the method. 
To demonstrate the efficiency of the extended loop al- 
gorithm, we apply the method to two different models, 
the Heisenberg antiferromagnetic model with easy-axis 
anisotropy and the Heisenberg spin ice model. For the 
former, we discuss the possibility of order-from-disordcr 
by thermal fluctuations. For the latter, we map out the 
phase diagram which includes a spin-ice like liquid state 
as well as a ferromagnetically-ordcred ice-rule state. 

This paper is organized as follows. In Sec. II, after 
reviewing the loop algorithm for Ising spin systems, we 
extend it to Heisenberg models with easy-axis anisotropy. 
In Sec. Ill and IV, we apply the extended loop algo- 
rithm to the pyrochlore antiferromagnet with the z-axis 
anisotropy and the spin-ice type fcrromagnct with (111) 
anisotropy, respectively. Summary is given in Sec. V. 



II. MONTE CARLO ALGORITHM 
A. Loop algorithm for Ising spin systems 

Before considering an extension of the loop algorithm 
to Heisenberg spin systems, here we briefly review the 
loop algorithm for Ising spin models with ice-rule type 
degenerac y 14 ' 17 . To generalize the following discussion, 
we assign black and white to two degrees of freedom of 
Ising spins in an appropriate manner. For example, for 
the spin ice model^, black and white represent inward 
and outward spins in terms of type-A tetrahedra, respec- 
tively, as shown in Fig. [TJ For the antiferromagnetic Ising 
model£, black and white simply correspond to up and 
down spins, respectively. Then the loop flip consists of 
two steps; first, we identify a closed loop which consists of 
alternating alignment of black and white sites, and then 
we try to flip all Ising spins on the loop. 

When all tetrahedra satisfy the ice rule as in the ice- 
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rule ground states, it is trivial to construct a loop of 
alternating black and white sites. At finite T, however, 
thermal fluctuations induce "defect tetrahedra" in which 
the "two-in two-out" ice-rule condition is violated and 
"three-in(out) one-out(in)" or "four-in(out)" configura- 
tion is realized (see Fig. [2]) . To maintain detailed balance 
in the loop flip, a loop must be contructed so that it does 
not involve defect tetrahedra on it. 

Two methods were proposed for the loop construction; 
one is called the long loop algorith m 15 i 17 and the other 
is the short loop algorithm^—. In the present study, we 
focus on the short loop algorithm. The procedure is the 
following (Fig. [2]): 

1. First, we randomly choose a tetrahedron which sat- 
isfies the ice-rule local constraint, namely, a tetra- 
hedron with two black and two white sites. 

2. We move to one of its four neighboring tetrahedra 
which also satisfies the ice rule, and mark the site 
shared by the two tetrahedra as the first site (the 
origin of path). If all neighboring tetrahedra arc 
defect, we abort the attempt to construct a loop 
and go back to the step I. 

3. Then we move to a tetrahedron which satisfies the 
ice rule out of three neighboring tetrahedra (except 
for the tetrahedra visited just before) , and mark the 
shared site as the second site. We repeat the pro- 
cedure and extend the path of marked sites which 
consists of an alternation of black and white sites. 
As in the step 2, once all three tetrahedra are de- 
fect, we abort the attempt and go back to the step 
1. 

4. If one of three neighboring tetrahedra has already 
been visited, we move to it and make a closed loop 
by deleting the dangling tail of the path. 

The procedure is slightly modified in the steps 3 and 4 
from the original short loop algorithm^—. In the orig- 
inal version, a move to neighboring tetrahedra is always 
completely random (except for the one visited just be- 
fore); it does not avoid defect tetrahedra (once a defect 
tetrahedron is chosen, the attempt is aborted) , and does 
not choose a previously-visited tetrahedron selectively. 
The above modifications enhance the efficiency by re- 
ducing the possibility to fail a loop construction, with 
satisfying the detailed balance in MC calculations. 

After the construction of a closed loop, all the col- 
ors on the loop are reversed simultaneously. This corre- 
sponds to a flip of all Ising spins on the loop: S, — > — Si. 
When the system has nearest-neighbor interactions only, 
the loop flip does not change the total energy; in other 
words, the flip is always accepted (rejection free update) 
in the MC sampling. When there are residual interac- 
tions such as farther-neighbor interactions, the loop flip 
is accepted with the probability depending on the total 
energy change by the standard Metropolis algorithm. 




FIG. 2. (color online). Schematic picture for a loop construc- 
tion by tracing a path through ice-rule tetrahedra. The path 
is made of alternating black and white sites. For simplicity, 
the figure shows a (lit) kagome layer with connected tetra- 
hedra. The path is denoted by a dashed line, and its lower 
part represents an example of a closed loop. 



In the long loop algorith m 15 ' 17 , a closed loop is formed 
only when the path returns to the initial site. This 
can generate a longer loop and its flip causes a bigger 
change of configurations. However, in general, it takes 
more CPU time to construct a longer loop. Moreover, a 
flip of a longer loop leads to a larger energy change and 
a lower acceptance rate. (For nearest-neighbor models, 
there is no energy cost and the long loop algorithm can 
be efficient.) Therefore, the short loop algorithm is more 
efficient than the long loop algorithm in general. 

At finite T, it is clear that the loop flip update does 
not satisfy ergodicity because it can change neither the 
spin configurations in defect tetrahedra nor the number 
of defect tetrahedra. It is, therefore, necessary to use 
the loop flip together with another update for retaining 
the ergodicity. This is easily achieved by introducing the 
standard single-spin flip in MC samplings. 



B. Extension of the loop algorithm to anisotropic 
Heisenberg spin systems 

In this section, we extend the loop algorithm to Heisen- 
berg models with easy-axis anisotropy. We start with a 
Hamiltonian of a general form: 

H = - JJ2 Si ■ Sj - Dj ]T (§ t ■ d,) 2 , (1) 
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where Si denotes a classical Heisenbcrg spin at site i on 
the pyrochlore lattice (we take |S,| = 1) and D\ (> 0) 
is the single-ion easy-axis anisotropy. The unit vector 
di defines the easy axis on site i. In the case of anti- 
ferromagnets with easy-axis anisotropy along the z axis, 
which are extensions of the Ising antiferromagneti, we 
set di = (0,0, 1) for all the sites. While, in the case of 
ferromagnets with the local (111) anisotropy, which are 
natural extensions of the spin ice models, we set a, to 
the direction connecting the centers of neighboring tetra- 
hedra from type B to A. 

To define a loop for the model given by Eq. (TTJ) , we as- 
sign black and white colors to sites at which Si-di > and 
Si ■ di < 0, respectively. Based on this definition, we can 
construct a closed loop with alternating black and white 
sites by following the short loop algorithm in Sec. II A. 
Then, we try to reverse all colors on the loop. However, 
this loop flip procedure is not unique in the Heisenberg 
spin case. To choose an efficient method, careful consid- 
eration on the energy change is necessary as discussed 
below. 

A natural extension of the loop flip in the Ising case 
is to reverse all three Cartesian components of Si on the 
loop as Si — > —Si, which we call flip xyz, as illustrated 
in Fig. [3] This flip xyz changes the energy by 

AE = 2J E $ ' & 

iGloop j'^loop 

= 2j E $ • &n + &o 

iGloop j'^loop 

= AE l{ +AE ± , (2) 

where 

A£|| =2J ]T Y, ^-S lh (3) 

zGloop j ^ loop 

AE ± = 2 J ]T E ^' ' &L> ( 4 ) 

igloop j'^loop 

and Si 1 1 (Sjj_) is the component of Si parallel (perpen- 
dicular) to di (see Fig. [3]). For models with nearest- 
neighbor interactions which retain ice-rule degeneracy in 
the ground state, one might expect that this update is 
always accepted in the limit of T — > 0; this is naively ex- 
pected since thermal fluctuations vanish and all the ice- 
rule configurations with spins parallel to the easy axes 
become energetically degenerate. However, this is not 
the case. At low T, spins fluctuate around the easy axes 
by angles of the order of T, and hence, the energy change 
by the flip is estimated as 

A£,| ex Y % « r2 > (5) 

iGloop 

AE ± cx Y °i K T > ( 6 ) 

iGloop 

where 8i denotes a deviation angle of spin % from 
the easy axis. Because the acceptance rate is 



given by min{l, exp(— AE/T)}, thermal fluctuations 
are irrelevant for the flip of Si || in the sense that 
lim^-vo cxp(— AEn /T) = 1. On the other hand, ther- 
mal fluctuations are relevant for the flip of Si± since 
liniT->.o ex p(— AE±/T) < 1. Therefore, flip xyz does not 
become rejection free in the limit of T — > even for 
nearest-neighbor models with ice-rule degeneracy. 

This consideration suggests a better way of a global 
update of spin configurations on the loop. That is a flip 
of only parallel components S^y as Si — > Si — 2(S, ■ di)di. 
This flip parallel changes the energy only by AE\\ cx T 2 
at low T, and hence, is expected to become rejection free 
as T — > for nearest-neighbor models with ice-rule de- 
generacy. The efficiency of these two updates, flip xyz 
and flip parallel, will be compared in numerical simula- 
tions in the following sections. 

As mentioned above, to retain the ergodicity, we use 
the loop flip together with the single-spin flip. One MC 
step consists of single-spin flips, followed by loop flips 
with either flip xyz or flip parallel. In the single-spin 
flips, we randomly choose a new spin state on the unit 
sphere for each spin following a procedure proposed by 
Marsaglia32. The loop flips are repeated until the number 
of tetrahedra visited exceeds the number of lattice sites. 
In our implementation, the procedure of loop flips takes 
CPU time comparable to that of a sweep of the lattice 
sites by single-spin flips. 

In the single-spin flips, we note that the completely 
random choice of a new spin direction leads to a low 
acceptance rate at low T, and a high acceptance rate 
is retained when restricting the new spin state within a 
small angle 8 around the original spin direction. How- 
ever, the high acceptance rate does not directly mean 
the high efficiency for the spin-ice type models, because 
such single-spin flips cannot retain the ergodicity at low 
T: small local fluctuations around the easy axes hardly 
lead to a global update between different spin-ice states 
at low T <C J- On the other hand, the random single- 
spin flips give a physically-important measure: its steep 
suppression signals formation of the spin-ice type man- 
ifold. This will be demonstrated for the two models in 
Sec. Ill and IV. 



III. APPLICATION TO PYROCHLORE 
HEISENBERG ANTIFERROMAGNETS WITH 
EASY-AXIS ANISOTROPY 

In this section, we apply the extended loop algorithm 
to the Heisenberg antiferromagnetic model with easy-axis 
anisotropy on the pyrochlore lattice. After introducing 
the model in Sec. Ill A, we demonstrate the efficiency 
of loop flips in MC simulations in Sec. Ill B. In Sec. 
Ill C, we discuss the possibility of order- from-disordcr 
phenomenon in comparison with related models. 
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(b) flip parallel ^ 



FIG. 3. (color online). Two different ways to flip black and 
white for the spin Si at site i: (a) flip xyz and (b) flip parallel. 
See the text for details. 



Model 



We consider an antiferromagnetic Hcisenberg model 
with easy-axis anisotropy along the z axis on the py- 
rochlore lattice. The model is given by taking <Jj = 
(0,0, 1) for all sites in the model ([1]). For simplicity, we 
take account of nearest-neighbor interactions only. The 
Hamiltonian is given by 



(7) 



(ij) 



Here we consider the antiferromagnetic case J < 0, and 
set an energy scale as |J| = 1, i.e., J = — 1. In the 
following MC calculations, we denote the linear dimen- 
sion of the system measured in the cubic unit cell by L. 
Namely, the total number of spins in the system N s is 
given by 16L 3 . Hereafter, we employ periodic boundary 
conditions. 

In the limit of Dj — > oo, this model reduces to the 
antiferromagnetic Ising model on the pyrochlore lattice. 
The ground state of the Ising model has the spin-ice 
type degeneracy, i.e., all spin configurations with two- up 
and two-down spins on every tetrahedron are energeti- 
cally degenerate^. It is known that the system does not 
show any phase transition at finite T in this case^i. For 
< D\ < oo, the situation at T = does not change; the 
ground state has the same macroscopic degeneracy. An 
interesting question is whether the degeneracy is lifted 
by thermal fluctuations in the Heisenberg spin model de- 
fined by Eq. ([TJ) . In other words, the question is whether 
an order-from-disorder phenomenon takes place in this 
model. To answer this question, it is necessary to inves- 
tigate the low-T properties, much lower than |J| where 
the ice-rule manifold is gradually formed. 



B. Demonstration of loop flip 

We investigate finite-T properties of the model ([7]) by 
MC simulation with the loop algorithm extended in Sec. 
II B. At low T compared to |J|, spin configurations are 
gradually enforced to satisfy the ice rule, and the accep- 
tance rate of the single-spin flip, P s ingic, is suppressed. 
This is demonstrated in Fig. @|a) for Di = 5.0: P S i n gic 
rapidly decreases at T ~ \J\ and vanishes almost linearly 
in T in the limit of T — > 0. On the contrary, the ac- 
ceptance rate of loop flips increases at low T. As shown 
in Fig. 0|a), the probability that a closed loop is suc- 
cessfully formed, Poop, steeply increases below T ~ |J|, 
indicating that almost all tetrahedra start to follow the 
ice rule below this temperature. At the same time, the 
acceptance rate of flips of a formed loop gradually in- 
creases at low T < | J\ and remains finite; here, P xyz and 
^parallel are the rate for flip xyz and flip parallel, respec- 
tively. The total acceptance rate of the loop flip is given 
by the product as Pi oop x P xyz or P 00 p x P par aiiei, and it 
sharply increases at T < \ J\, compensating the decrease 

of Psingic ■ 

As clearly indicated in Fig.^Ja), the acceptance rate of 
flip parallel is always larger than that of flip xyz. In par- 
ticular, Pparaiid approaches 1 (rejection free) as T — >• 0, 
whereas P xyz goes to a smaller value ~ 0.5. The reduc- 
tion of P xyz becomes larger for smaller anisotropy D\. 
This is demonstrated at T = 0.1 in Fig. BJb); P xyz de- 
creases almost exponentially with 1/P>i and decreases to 
about 0.01 at D\ = 0.5. On the other hand, P pa raiici is 
almost independent of D\ and remains rejection free as 
T — > in the wide range of Dj. Note that the single-spin 
flip does not work efficiently even for weak anisotropy, 
e.g., for Di = 0.5, P sing i ~ 0.03 at T = 0.1. There- 
fore, flip parallel retains higher efficiency than flip xyz 
and compensates the low efficiency of the single-spin flip 
over a wide range of D\. These behaviors are consistent 
with the argument in Sec. II B, and demonstrate the 
advantage of flip parallel at low T. 

To further demonstrate efficiency of the loop flips, we 
calculate the autocorrelation function of spin configura- 
tions. Here the autocorrelation function is defined for an 
interval of n MC steps in the form 



A(n) = — 



2J & (n ) • Si(n + n) 



(8) 



where N s is the number of lattice sites. The results at 
T = 1.0,0.5,0.05 are shown in Fig. [SJ We average the 
data over independent 10 4 runs after no = 1 x 10 4 ther- 
malization. The decay of the autocorrelation function for 
the single-spin flip dynamics rapidly becomes slower at 
low T: The autocorrelation remains almost 1 for n < 100 
at the lowest T = 0.05. This clearly indicates the freez- 
ing of single-spin flip at low T. In contrast, the loop-flip 
dynamics exhibits no signature of freezing down to the 
lowest T in the present calculations. In particular, even 
at T = 0.05, the autocorrelation vanishes rapidly for flip 
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FIG. 4. (color online), (a) Temperature dependence of the 
acceptance rate of the single-spin flip (P S mgie), the probability 
of formation of closed loops (Pioop), the acceptance rates of 
flip of a formed loop by flip xyz {P xyz ) and by flip parallel 
(-Pparaiiei). The data are calculated for the model ([7]) at D\ — 
5.0. The data for the system sizes L — 2 and L = 4 are 
denoted by crosses and filled squares, respectively, (b) D\ 
dependence of P pa raliei, Pxyz, Poop, and P s i ng ic at T = 0.1 for 
L = 2. 



parallel; it becomes smaller than 0.01 for n > 3, indi- 
cating the autocorrelation time is 0(1). On the other 
hand, as clearly indicated in the inset of Fig. [SJc), the 
loop-flip dynamics with flip xyz exhibits a slower residual 
relaxation for n > 10. This slow relaxation may be due 
to enhanced spin fluctuations around thermally-induced 
"defect tetrahedra" , which lower the efficiency of flip xyz. 
The flip parallel does not severely suffer from such fluc- 
tuations and keeps the high efficiency down to the lowest 
T. 



C. Absence of order from disorder 

As mentioned in Sec. Ill A, in the Ising limit D\ — > oo, 
the model does not show any phase transition, and the 
ground state has continuous macroscopic degeneracyi^i. 
On the other hand, in the Heisenberg limit, i.e., D\ = 0, 
the model shows no transition down to T = and the 
ground state has a continuous degeneracy^—. However, it 
was pointed out that, for the planar XY type anisotropy 



(a) 



(b) 



(c) 



1 



r f= \ flip parallel 
flip xyz 
single-spin flip 




80 100 



n 



FIG. 5. (color online). Autocorrelation functions for spin 
configurations [Eq. ©] calculated at (a) T = 1.0, (b) T = 0.5, 
and (c) T = 0.05. The data are at Dj = 5.0 for L = 2. The 
offset ^(oo) was obtained by averaging A(n) in the range of 
5000 < n < 10000. The inset of (c) shows the same data as 
in (c) in a wider range of n. 



with D\ < 0, although the system suffers from a continu- 
ous degeneracy in the ground state, thermal fluctuations 
select a subset from the continuous manifold and induces 
a first-order transition to a conventional Neel order—. 
Furthermore, on the 2D kagome lattice (a (111) plane in 
the pyrochlore lattice), XX Z models with Ising- type ex- 
change anisotropy also exhibit a thermally-driven phase 
transition by selecting a subset from an ice-rule type dis- 
crete manifold of the ground stated. These results are 
examples of the so-called order- from-disorder phenomena 
appearing in the systems with anisotropy. 

Motivated by these previous studies, we here exam- 
ine the possibility of order-from-disordcr in the present 
model with Ising type anisotropy < D\ < 00. We 
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calculate the specific heat C and the uniform magnetic 
susceptibility xo by MC simulation with the single-spin 
flip and the loop flip of flip parallel. The specific heat C 
and the uniform magnetic susceptibility \o are calculated 
by 



C 



1 (H 2 ) - (H)' 



T 



and 



Xo 



1 (A/ 2 ) 
3ivT T 



(9) 



(10) 



respectively, where (• • ■) denotes a thermal average. Here 
the square of total magnetization M 2 is given by M 2 = 
Yl,fj,=x y z (Si Si) 2 ■ Numbers of MC steps for thermal- 
ization, N th _, and for sampling, iV samp , are (AT th , iV samp ) = 
(1 x 10 4 , 1 x 10 5 ) for L = 2 and L = 4, (5 x 10 4 , 5 x 10 5 ) 
for L = 5, respectively The data are averaged over four 
independent MC runs to estimate statistical errors by 
variance of average values in the runs. 

Figure [B] shows the results at Dj = 5.0. Although we 
find a broad peak in the specific heat C at T ~ 0.85 
corresponding to cooperative formation of ice-rule tetra- 
hedra [see also Pi 00 p in Fig. 0], there is no signature 
of a phase transition in the specific heat and the mag- 
netic susceptibility down to the calculated lowest tem- 
perature T = 0.02. The data do not show any sin- 
gularity and collapse onto each other among different 
sizes for L < 5. This indicates the absence of ordcr- 
from-disorder in the present model with Ising anisotropy 
D\ > 0. We confirmed that the situation is unchanged for 
several other values of positive D\ . The absence of order- 
from-disorder is seemingly consistent with a Maxwellian 
counting argument as follows^. The pyrochlore Heisen- 
berg model is more underconstrained than the kagome 
one which exhibits ordcr-from-disorder. In fact, the 
Heisenberg model on the pyrochlore lattice is a marginal 
model in terms of ordcr-from-disorder; namely, a simple 
Maxwellian counting alone cannot conclude whether the 
order- from-disorder occurs or not. Our numerical results 
show that it does not occur in the present model with 
easy-axis anisotropy, in contrast to the case of easy-plane 
anisotropy. This contrasting behavior depending on the 
sign of D\ appears to support that the pyrochlore case is 
marginal. 



IV. APPLICATION TO HEISENBERG SPIN ICE 
MODEL 

In this section, we apply the extended loop algorithm 
to the Heisenberg spin ice model. After introducing the 
model in Sec. IV A, we demonstrate the efficiency of 
loop flips in MC simulation in Sec. IV B. In Sec. IV C, 
we clarify the phase diagram, and discuss the nature of 
crossover and phase transition among paramagnet, spin- 
ice liquid, and fcrromagnctically-ordered spin-ice state. 
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FIG. 6. (color online). Temperature dependences of (a) the 
specific heat C and (b) the uniform magnetic susceptibility 
Xo for the model at D\ — 5.0. 



A. Model 

We consider a Heisenberg ferromagnet with the local 
(111) easy-axis anisotropy (Heisenberg spin ice model), 
whose Hamiltonian is given by 



H 



(11) 



Here Si denotes a classical Heisenberg spin at site i and 
D\ (> 0) is the single- ion easy-axis anisotropy. The unit 
vector di defines the local (111) easy axis on site i, which 
is along the direction connecting the center of two tetra- 
hedra sharing the site i from type B to A (see Fig. [IJ. 
For simplicity, we take account of nearest-neighbor in- 
teractions only. We consider the ferromagnetic exchange 
J > 0, and take the energy unit as J = 1. As in the 
previous section, the system size is N s = 16L 3 spins in 
the following calculations. 

The Ising counterpart of this model (D\ = oo) is the 
nearest-neighbor spin-ice model, whose ground state re- 
tains the ice-rule degeneracy. The model does not ex- 
hibit any phase transition down to T = 0. However, it 
shows a crossover at T* ~ O(J) related with cooperative 
formation of the ice-rule degenerate manifold, which is 
signaled by a broad peak in the specific heat^^. When 
< Di < oo, the system tends to gain the exchange 
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energy by canting spins from the easy axes, as schemat- 
ically shown in Fig. [TJa). Through this canting, the net 
moments in each tetrahedron associated with the two-in 
two-out spin configuration are aligned among tetrahedra 
to minimize the exchange energy; therefore, we expect a 
canted ferromagnetic ground state for < D\ < oo, as 
shown in Fig. [TJb). 

In fact, a recent MC study showed that the model with 
a positive D\ exhibits a finite-T phase transition^: The 
system exhibits a transition from the high-T paramagnet 
to the low-T ferromagnctically-ordcred ice-rule state, and 
the transition temperature T c increases as D\ decreases. 
However, the obtained phase diagram is limited to a re- 
gion with relatively- weak anisotropy, D\ < 25. This is 
because the MC study was done by using the single- 
spin flip algorithm which suffers from a severe freezing 
in strongly anisotropic cases. In the calculated region of 
D\, the transition temperature T c is always larger than 
T* , and the ordering takes place before forming the ice- 
rule manifold at T* . 

Hence, there remains an interesting question: What 
happens when the ice-rule manifold is formed before the 
ordering? In other words, how does the formation of ice- 
rule manifold affect the phase transition? In the vicinity 
of the Ising limit D\ ^> J, we expect that T c becomes 
very small and even lower than T* . Hence, the system 
provides a chance to investigate the relation between a 
crossover from high-T paramagnet (gas) to a cooperative 
ice-rule state (liquid) [Fig. [7Jc)] and the phase transition 
from the spin-ice liquid to the ferromagnetic ordering of 
the ice- ruled tetrahedra (solid). This is a gas-liquid-solid 
transition in terms of spins. To elucidate their nature, 
it is necessary to properly sample the ice-rule manifold 
without dynamical freezing, which will be a good target 
of the extended loop algorithm. 



B. Demonstration of loop flip 

We demonstrate the efficiency of the loop algorithm 
for the strong easy-axis anisotropy, i.e., D\ ^> J. Fig- 
ure |UJa) shows the comparison of the acceptance rates 
at D\ = 50.0 as an example. Because of the large 
Dj, the acceptance rate of the single-spin flip, P S i ng i e , 
is strongly suppressed and becomes less than 1% already 
at T ~ J (cf. Fig. 01); with further decreasing T, P S mgle is 
steeply reduced because of gradual formation of the ice- 
rule manifold. On the contrary, the probability of loop 
formation, Pi 00 p, rapidly increases below T ~ J. At the 
same time, the acceptance rates of flip parallel and flip 
xyz, -Pparaiiei and P X yz, remain much larger compared to 
-fsingie- Hence, the loop flips are effective down to low 
T and compensate the freezing of the single-spin-flip dy- 
namics. It is also noted that Pparaiiei is always larger 
than P X y Z , being consistent with the argument in Sec. II 
B. In contrast to the model in Sec. Ill, however, P pa raiiei 
and P xy z do not approach a finite value as T — > but 
exhibit a sharp drop at T ~ 0.16. This is due to the fer- 




(b) Ferromagnetically-ordered state 




(c) Spin-ice liquid 




FIG. 7. (color online), (a) Schematic picture for a ground 
state of the Heisenberg spin ice model given by Eq. {TT} on 
an isolated tetrahedron. 8 C is a canting angle. The big ar- 
row in the center of tetrahedron denotes a net moment by 
forming the canted two-in two-out spin configuration, (b) 
A ferromagnetically-ordered ice-rule state with aligning net 
magnetizations along the z axis, (c) Schematic picture for a 
spin-ice liquid, in which each tetrahedron is in the two-in two- 
out state but the induced net moments are disordered among 
tetrahedra. 



romagnetic transition mentioned above. We will analyze 
the nature of the transition in more details in the next 
subsection. 

To further demonstrate the efficiency of the loop flips 
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FIG. 8. (color online), (a) Temperature dependence of the 
acceptance rates of different updates. The notations are the 
same as in Fig. [4] (b) MC dynamics of the squared uniform 
magnetization m 2 in the MC runs with and without flip par- 
allel at T = 0.05. The data are for L — 2. Two data for 
each case show the results starting from different initial spin 
configurations. 



at low T, we calculate MC dynamics of the squared mag- 
netization m 2 for L = 2 at T = 0.05 (below the transition 
temperature T c ). The squared magnetization m 2 is de- 
fined by 



2 _ 



(12) 



As shown in Fig. EJb), when we employ the single-spin 
flip only, MC dynamics shows a severe freezing; m 2 does 
not converge to an expected value ~ 0.37 (the value will 
be discussed later) even after 4 x 10 5 MC steps, and 
moreover, the data are frozen after ~ 10 4 MC steps at 
some different values depending on initial spin config- 
urations. On the other hand, MC dynamics is greatly 
accelerated by the loop flip and m 2 converges to the ex- 
pected thermal-equilibrium value after 5 x 10 4 MC steps. 
These clearly show the advantage of the loop algorithm 
in investigating the low-T properties of the Hciscnberg 
spin ice model for D\ J. 



C. Spin-ice liquid-to-solid transition 

We investigate thermodynamic properties of the 
Heisenbcrg spin ice model using the loop flip of flip par- 
allel and the single-spin flip. To accelerate the MC dy- 
namics further, particularly at very low T where P pa raiiei 
is suppressed, we adopt the overrelaxation update^ and 
the exchange MC method 2 ^. In the overrelaxation up- 
date, the energy change arising from the second term 
in Eq. (fTTj) is treated by the standard Metropolis algo- 
rithm. Numbers of MC steps for thermalization, iV t h, 
and for sampling, iV samp , are (iV th; iV samp ) = (5 x 10 5 , 5 x 
10 5 ), (5 x 10 5 , 5 x 10 5 ), (2 x 10 6 , 2 x 10 6 ) for L = 2,3, 4, 
respectively. The data are averaged over four indepen- 
dent MC runs to estimate statistical errors by variance 
of average values in the runs. 

First, we show T dependences of the squared mag- 
netization m 2 and the specific heat C at D\ = 50.0 in 
Fig. [SJ As shown in Fig. Ela), rn? exhibits a steep rise 
at T ~ 0.16. This rise becomes steeper with increas- 
ing system size and almost discontinuous in the largest 
size L = 4. The saturation value is m? ~ 0.37 which is 
expected for the canted ferromagnetic ordered state dis- 
cussed in Sec. IV A (the value will be discussed below). 
At the same temperature, C exhibits a sharp peak, which 
also becomes almost discontinuous at L = 4. These in- 
dicate that the system exhibits a first-order transition to 
the ferromagnetic state with a canted ice-rule spin con- 
figuration [Fig. Efb)]. These sharp singularities allow us 
to determine the transition temperature T c with enough 
accuracy for the present purpose. 

Another anomaly is found as a broad peak of C at T* ~ 
0.25 as shown in the inset of Fig. [HJb). Since the broad 
peak does not show any significant size dependence and 
since a similar peak is also seen in the Ising counterpart, 
this is a crossover from paramagnct to a cooperative ice- 
rule state; below T*, all tetrahedra tend to satisfy the 
ice-rule spin configurations but the system still remains 
paramagnetic for T > T c . In fact, T* coincides with the 
sharp rise of Pi oop in Fig. HJa). We call this correlated 
state for T c < T < T* spin-ice liquid, whose schematic 
picture is shown in Fig.JTJc). 

Now we examine D\ dependence of the phase tran- 
sition and the crossover. Figure [TU] shows the cal- 
culated results of m 2 and C for L = 4 at A = 
100.0, 50.0, 33.3, 25.0, 16.6. Numbers of MC steps for 
thermalization, N t h, and for sampling, iV samp , are 
(^th^samp) = (1 x 10 6 ,1 x 10 6 ) for Dj = 100.0,25.0, 
(4 x 10 6 ,4 x 10 6 ) for D x = 33.3, and (6 x 10 4 ,6 x 10 4 ) 
for £>i = 16.6. As shown in Fig. [TU1 T c , signaled by a 
sharp rise of m 2 and a singular peak in C, increases as 
D\ decreases. On the other hand, T*, at which C shows 
a broad peak, does not strongly depend on D\. As a 
consequence, for 1/Dj > 0.04, T c becomes higher than 
T*; the spin- ice liquid state is restricted to T c < T < T* 
for \/D\ < 0.04. Note that the previous MC study was 
restricted to the region of l/Dj > 0.04 where T c > T* 
and the single-spin flip does not show severe freezing by 
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FIG. 9. (color online). Temperature dependences of (a) the 
squared uniform magnetization m 2 and (b) the specific heat 
C at Di — 50.0. A broad peak in the specific heat is denoted 
by an arrow. The inset shows an enlarged view of the main 
panel. 



the formation of the ice-rule manifold at T*. 

At the lowest T, m 2 approaches a constant whose value 
is dependent on Dj. By considering the ground-state 
energy following the discussion in Ref. [27l the saturation 
value can be calculated as 



m 2 {T = 0) = - (V2sin(9 c + cos0, 



2 

c I ) 



where 8 C is the optimal canting angle 

I 4 



T 1 

u c = arctan ■ 

4 2 16 V J 



(13) 



(14) 



Note that m 2 = in the Ising limit (1/Dj = 0), while 
rn 2 —> 1/3 in the limit of 1/D\ 0: The anisotropy 
D\ is a singular perturbation to the macroscopically- 
degenerate ground-state in the Ising case. As shown in 
Fig. rTUT c). to 2 at the lowest T in our MC simulation scales 



well with Eq. (fT3"]) : This demonstrates the efficiency of 
our algorithm down to the lowest T even in the large D\ 
region where T c < T* . 

We summarize the results in the phase diagram in 
Fig. QTJ As mentioned above, T c grows as 1/Di in- 
creases; in particular, it almost linearly increases in the 
small 1/Di region. On the other hand, T* remains al- 
most constant irrespective of D\. Consequently, in the 




0.06 



FIG. 10. (color online). Temperature dependences of (a) the 
squared uniform magnetization m 2 and (b) the specific heat 
C. (c) D\ dependence of m 2 at the lowest temperature in 
(a); the dotted line denotes the ground-state magnetization 
calculated from Eq. (|13[) . The system size is L — 4. 



region of 1/D\ < 0.04, there is a successive crossover and 
phase transition, i.e., a crossover from the high-T para- 
magnet (gas) to the intermediate-T spin-ice liquid at T* , 
and a transition from the spin-ice liquid to the low-T 
ferromagnetically-ordcrcd spin- ice solid at T c . In con- 
trast, for 1/Dj > 0.04, there is a single transition from 
the high-T paramagnet to the low-T spin-ice solid, being 
consistent with the previous result. The transition at 
T c is of first order, while the discontinuity appears to be- 
come weaker as l/D\ increases^ (see Fig. ITU)) . The phase 
diagram illuminates the rich physics in the Heisenberg 
spin ice model, including the interesting spin-ice liquid- 
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FIG. 11. (color online). Di-T phase diagram of the Heisen- 
berg spin ice model given by Eq. (|11|> . The transition tem- 
perature T c and the crossover temperature T* are estimated 
from a sharp peak and a broad one in the specific heat for 
L = 4, respectively. 



to-solid phase transition, which has not been revealed in 
the previous study. 

A similar liquid-to-solid transition was recently dis- 
cussed for the bilinear-biquadratic Heisenberg model on 
the pyrochlore lattice in applied magnetic fieldii. In this 
model, a strong biquadratic interaction under the influ- 
ence of magnetic field enforces the ice-rule like local con- 
straint on spin configurations. Farther-neighbor interac- 
tions lift the ice-rule like degeneracy and induce a phase 
transition to a long-range ordered state. The phase dia- 
gram is quite similar, and the transition is of first order 
also in this case. 



parallel becomes rejection free as T — > but flip xyz not 
for models in which the ground state has macroscopic 
ice-rule degeneracy. 

We have demonstrated the efficiency of the loop flips by 
performing MC simulations for two typical models, the 
Heisenberg antiferromagnet with easy-axis anisotropy 
along the z axis and the Heisenberg spin ice model. By 
using the extended loop algorithm, we have investigated 
low-T properties of the two models, which are hard to ac- 
cess by the standard single-spin flip alone. For the former 
model, we have critically checked the absence of order- 
from-disorder phenomenon. For the latter model, we 
have successfully obtained the rich phase diagram involv- 
ing the gas-liquid-solid like transition among the param- 
agnet, spin-ice liquid, and ferromagnetically-ordered ice- 
rule state. 

Before closing this paper, we make a brief remark on 
the application of the extended loop algorithm to the 
bilinear-biquadratic Heisenberg model on the pyrochlore 
lattice. It was recently pointed out that ferromagnetic bi- 
quadratic interactions lead to formation of ice-rule type 
manifold and farther-neighbor bilinear interactions per- 
turb this manifold to select an ordered stated. In par- 
ticular, under applied magnetic field, the system shows a 
gas-liquid-solid like transition, similar to the Heisenberg 
spin ice model discussed in this paper. The low-T prop- 
erties have not yet been fully clarified, because of the 
freezing of single-spin flips. We expect that our extended 
loop algorithm works efficiently also in this model. Such 
extension will be reported elsewhere. 
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V. SUMMARY 
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